%{
Calculate the annual and seasonal mean of the variable for each year
%}

%%% define period
period = 'annual';  %   'annual'    'MAM'  'JJA'  'SON'   'DJF'

disp(period)

%%% define the file name of all models 
model = { 'albm', 'lake', 'simstrat-uog'};
gcm = { 'gfdl-esm2m', 'hadgem2-es', 'ipsl-cm5a-lr', 'miroc5' };
rcp = { 'historical', 'picontrol', 'rcp26', 'rcp60', 'rcp85' };


tic

for flag_model = [1:4]      %   #1 = albm      #2 = lake      #3 = simstrat-uog    

%%% define target timestep by define start and end time (annual/seasonally)
if strcmp( period,'annual' )
    
    time_serises_sta = { datetime([1861:2005], 1, 1), datetime([1661:1860], 1, 1), ...
                    datetime([2006:2099], 1, 1), datetime([2006:2099], 1, 1), datetime([2006:2099], 1, 1) };   % for 4 rcp
    time_serises_end = cellfun( @(x) x + calmonths(12) - days(1), time_serises_sta, 'UniformOutput', false );
    
elseif strcmp( period,'MAM' )

    time_serises_sta = { datetime([1861:2005], 3, 1), datetime([1661:1860], 3, 1), ...
                datetime([2006:2099], 3, 1), datetime([2006:2099], 3, 1), datetime([2006:2099], 3, 1) };   % for 4 rcp
    time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false );
    
elseif strcmp( period,'JJA' )
    
    time_serises_sta = { datetime([1861:2005], 6, 1), datetime([1661:1860], 6, 1), ...
                    datetime([2006:2099], 6, 1), datetime([2006:2099], 6, 1), datetime([2006:2099], 6, 1) };   % for 4 rcp
    time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false );
    
elseif strcmp( period,'SON' )
   
    time_serises_sta = { datetime([1861:2005], 9, 1), datetime([1661:1860], 9, 1), ...
                    datetime([2006:2099], 9, 1), datetime([2006:2099], 9, 1), datetime([2006:2099], 9, 1) };   % for 4 rcp
    time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false );
    
else
   
    time_serises_sta = { datetime([1861:2004], 12, 1), datetime([1661:1859], 12, 1), ...      % 今年的1、2月和去年的12月算今年的冬季，因此起始年无冬季数据
                    datetime([2006:2098], 12, 1), datetime([2006:2098], 12, 1), datetime([2006:2098], 12, 1) };   % for 4 rcp
    time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false ); 
    
end
    
%% Daily to annual or seasonaly 

    for flag_gcm = 1:4      %   #1 = gfdl-esm2m      #2 = hadgem2-es      #3 = ipsl-cm5a-lr      #4 = miroc5

        for flag_rcp = 1:5      %    #1 = historical    #2 = picontrol      #3 = rcp26      #4 = rcp60      #5 = rcp85
            
            %%% var. to save
            surftemp_mean = [];
            source = [];
            
            %%% change time serises for albm in picontrol
            if flag_model == 1 && strcmp( period,'annual' )
                time_serises_sta{2} = datetime([2006:2099], 1, 1);   % for 4 rcp
                time_serises_end = cellfun( @(x) x + calmonths(12) - days(1), time_serises_sta, 'UniformOutput', false );
            elseif flag_model == 1 && strcmp( period,'MAM' )
                time_serises_sta{2} = datetime([2006:2099], 3, 1);   % for 4 rcp
                time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false );
            elseif flag_model == 1 && strcmp( period,'JJA' )
                time_serises_sta{2} = datetime([2006:2099], 6, 1);   % for 4 rcp
                time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false );
            elseif flag_model == 1 && strcmp( period,'SON' )
                time_serises_sta{2} = datetime([2006:2099], 9, 1);   % for 4 rcp
                time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false );
            elseif flag_model == 1 && strcmp( period,'DJF' )
                time_serises_sta{2} = datetime([2006:2098], 12, 1);   % for 4 rcp
                time_serises_end = cellfun( @(x) x + calmonths(3) - days(1), time_serises_sta, 'UniformOutput', false ); 
            else
            end
            
            for flag_time = 1:flag_time_e
                    
                %%% file name for input(surftemp data)        
                loc_input_surftemp = [ '/home/ISIMIP/output/', model{flag_model}, ...
                            '/', gcm{flag_gcm}, '/',rcp{flag_rcp}, '/' ];
                file_input_surftemp = [ model{flag_model},   '_', gcm{flag_gcm}, '_', rcp{flag_rcp}, '_', 'surftemp_daily_', ...
                            num2str( year( time_serises_sta{ flag_rcp }(flag_time) ) ), '.mat' ];
                
                %%% load surftemp
                load( [loc_input_surftemp, file_input_surftemp], 'surftemp' ) % load surftemp, var name is 'surftemp'
                
                % define the start and end location
                sta_loc = days( time_serises_sta{ flag_rcp }(flag_time) - ...
                          datetime( year( time_serises_sta{ flag_rcp }(flag_time) ),1 ,1 ) ) + 1;
                end_loc = days( time_serises_end{ flag_rcp }(flag_time) - ...
                          datetime( year( time_serises_end{ flag_rcp }(flag_time) ),1 ,1 ) ) + 1;
                
                 % cut data 
                if strcmp( period,'DJF' )   % DJF requires cross-file reading（December for the first file plus January and February for the next file）  
                    
                    surftemp_cut = surftemp( :, :, sta_loc:end ) ;   
                    
                    % load the next year data for month 1 2
                    file_input_surftemp = [ model{flag_model},   '_', gcm{flag_gcm}, '_', rcp{flag_rcp}, '_', 'surftemp_daily_', ...
                                num2str( year( time_serises_sta{ flag_rcp }(flag_time) ) + 1 ), '.mat' ];

                    load( [loc_input_surftemp, file_input_surftemp], 'surftemp' ) % load surftemp, var name is 'surftemp'
                    
                    surftemp_cut = cat( 3, surftemp_cut, surftemp( :, :, 1:end_loc ) ) - 273.15;    % Convert to Celsius
                    
                else
                    surftemp_cut = surftemp( :, :, sta_loc:end_loc ) - 273.15;    % Convert to Celsius
                end
                
                %%% show the year and sta/end loc
                disp( [ model{flag_model},   '_', gcm{flag_gcm}, '_', rcp{flag_rcp}, '_', 'surftemp_daily_', ...
                            char( time_serises_sta{ flag_rcp }(flag_time),'yyyy' ), ' ', num2str(sta_loc), '-' num2str(end_loc) ] )
                
                %%% is lake
                notlake = isnan( sum( surftemp_cut, 3 ) ); % Identify non-lake points
                 
                %%% surftemp mean
                surftemp_mean_timestep = mean( surftemp_cut, 3, 'omitnan' );
                surftemp_mean_timestep( notlake ) = nan;
                
                %%% save to the period
                surftemp_mean(:, :, flag_time) = surftemp_mean_timestep;
                
                source ( flag_time, 2 ) = sta_loc;
                source ( flag_time, 3 ) = end_loc;
                
                if strcmp( period,'DJF' )   % DJF is counted as the year in which January and February are located, so +1 
                    source ( flag_time, 1 ) = year(time_serises_sta{ flag_rcp }(flag_time)) + 1;
                else
                    source ( flag_time, 1 ) = year(time_serises_sta{ flag_rcp }(flag_time));
                end
                
            end
            
            %%% DJF is counted as the year in which January and February are located
            if strcmp( period,'DJF' )   % No DJF data for the first year of the time series, set to nan
                inter_nan = ones( size(surftemp_mean,1), size(surftemp_mean,2) );     
                inter_nan( inter_nan == 1 ) = nan;
                surftemp_mean = cat( 3, inter_nan, surftemp_mean );
                source = [ [nan, nan, nan]; source ];
            else
            end
            
            %%% save to mat file
            loc_output = [ '/home/work_file/', model{flag_model}, '/' ];
            file_output = [ period, '_', model{flag_model},   '_', gcm{flag_gcm}, '_', rcp{flag_rcp}, ...
                            '_LSWT', '.mat' ];
            save( [ loc_output, file_output ], 'surftemp_mean','source' )

        end

    end

end

toc





